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AXISYMMETRICAL  CREEP  BUCKLING  OF  CLAMPED 
SHALLOW  SPHERICAL  SHELLS* 

N.  C.  Huang 
Stanford  University 

Introduction 

The  problem  of  axlsymmetrlcal  buckling  of  clamped 

shallow  spherical  elastic  shells  under  uniform  external 

/  \  fl-41** 

pressure  (Figure  1)  has  been  studied  by  many  authors 

with  different  numerical  methods.  Consistent  results  for 

the  buckling  pressure  for  shells  with  a  large  range  of 

geometrical  parameters  were  obtained  in  their  work.  In 

this  paper  we  shall  consider  the  effect  of  viscoelastic 

properties  of  shells  on  the  axlsymmetrlcal  buckling.  For 

viscoelastic  shells  under  constant  uniform,  pressure,  the 

average  vertical  deflection  increases  with  time  as  shown 

in  Figure  2.  Whan  time  reaches  a  critical  value  t^^  , 

the  rate  of  increase  of  the  average  vertical  deflection 
becomes  very  high,  which  indicates  the  initiation  of 
snapping  of  the  shell  at  this  Instant.  After  snapping, 
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the  shell  continues  to  deform  as  shown  in  Figure  2.  Our 
interest  is  to  find  the  critical  time  t^^  for  different 

pressures  and  for  shells  with  different  geometrical  para¬ 
meters.  The  pressure  on  shells  is  considered  to  be 
applied  either  by  pressurized  air  in  a  chamber  of  large 
volume  or  by  an  incompressible  fluid.  For  the  latter 
loading  condition,  the  average  deflection  remains  constant 
and  the  applied  pressure  relaxes.  In  this  paper,  the 
relaxation  of  external  pressure  of  this  type  is  analyzed, 
as  well  as  prescribed  pressure  increases. 

Basic  Equations 

Consider  an  elastic  clamped  shallow  spherical  shell 
as  shown  in  Figure  1.  For  axi symmetrical  deformation, 
the  governing  differential  equations  and  boundary  con¬ 
ditions  are  found  in  [1]  as 
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where  6  =  -  j—  j  w  is  the  vertical  deflection  of  the 

shell,  f  is  a  stress  function,  and  E  and  v  are  the 

Young's  modulus  and  Poisson's  ratio  respectively.  Let 

and  be  two  constants  and  put 
o  o  ^ 
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Equations  (1-4)  can  then  be  written  as 

i2(x^e"  +  X  O'  -  e)  =  -  x^<t>  -  2px^  +  X  e  0 

X  <i>"  +  X  If'  -  D  =  (x^e  -  i  xe^) 


(5) 
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(6) 


0(A)  =  0 


(7) 


A<I>'(A)  -  V  <t(A)  =  0  (8) 

where  (  )  =  35^  (  )  •  The  material  properties  are  involved 
in  the  constants  and  Ig. 

In  our  viscoelastic  theory^  we  use  the  following 
assumptions:  (l)  the  strain  components  in  the  shell  are 
small  so  the  constitutive  equations  are  linear;  (2)  for 
mathematical  convenience,  the  bulk  modulus  of  the  shell,  K, 
is  assumed  to  be  independent  of  time;  (3)  the  Kirchhoff's 
hypotheses  hold  in  the  whole  process  of  deformation; 

(4)  the  shell  is  assumed  to  deform  quasi-statically,  i.e. 
the  effect  due  to  its  inertia  force  is  neglected. 

In  the  viscoelastic  problem,  the  Young's  modulus  E 
is  replaced  by  an  integral  operator  with  a  kernal  function 
E(t).  The  kernal  function  equivalent  to  Poisson's  ratio 
in  elastic  case,  v(t),  can  be  found  directly  from  the 
relation 


v(t) 


1  _  Elt) 

2  6K 
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We  use  the  notation  =  E(0)  and  =  v(0),  then 

q^  defined  above  is  the  classical  buckling  pressure  of 

a  complete  spherical  shell  with  the  sam.e  thickness  and 
radius  of  curvature.  The  constants  I^  and  I2  are 

replaced  by  the  integral  operators,  with  kernal  functions 
l^(t)  and  l2(t)  which  are  determined  by  the  following 
equations : 
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where  the  starred  notation  represents  the  Integral  operation 
as  shown  In  the  following  example: 


I  *  f 


t 


—00 


Thusj  equation  (11)  Is  a  Volterra  Integral  equation  for  I^. 
Examining  the  derivation  of  equations  (l-^)  In  [1],  we 
are  able  to  derive  the  governing  equations  for  the  defor¬ 
mation  of  viscoelastic  shells  under  the  above  assumptions 
by  applying  the  correspondence  principle  to  equations 
(5-8): 


I*  (x^0"+  x0'-e)  =  -  x^C)  -  2px3  +xe<l)  (12) 

xO"  +  Xi>'  -  ^  =  l^*(x^0  -  ^X0^)  (13) 

e{\)  =  0  (14) 

AO’(A)  -  v*4>(?v)  =  0  (15) 
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where  6,  $  and  p  are  time  dependent  and  (  ) '  =  ^  (  )  . 
Equations  (12)  and  (13)  are  two  nonlinear  Integro- 
dlfferentlal  equations  for  the  two  unknown  functions 
e(xjt)  and  $(x,t).  These  equations  will  be  solved  numerical¬ 
ly.  The  average  vertical  deflection  p  is  found  in  [1]  as 
A 


x2 
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Let  M  and  M, 

r  X/ 

moments  respectively. 
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From  the  stress  strain  relations^  we  have 

"■r-  V  -  I3  *  <f) 

"■t  ■  ^2  '  -  I3  ,  e'  (20) 

where  the  operator  is  defined  as 

I3  =  I2  *  V  (21) 

Hence,  and  can  be  obtained  from  9  by  numerical 
differentiation  and  integrations. 
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Numerical  Method 


The  viscoelastic  material  in  our  investigation  is 
polymethylmethacrylate.  The  nondimensional  relaxation 
modulus  is  found  in  [5]  as 

I^(t)  =  10  +  erf[0.3l(log^Qt  -  3.6)]}  (22) 

we  take  =  0.35  which  is  equivalent  to  K/E^  =  1.111. 
Let  t  =  t^.  We  shall  use  the  following  approximation  as 
shown  in  [6]  for  the  time  integration: 

I*f  =  l(0)f^  -  i  V  {(fj  ^ 

=  I  [f^  G(l,t^)  -  F(l,f^)]  (23) 

where  F(l,f^)  =  ‘  ^  (^n"^n-l)] 

(24) 

G(l,t^)  =  1(0)  +  (25) 

From  the  above  definition,  F(l,fj^)  is  linear  with 
respect  to  f^^,  i.e. 
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F(I>  +  Fd^S^) 

F(l,kfj^)  =  kF(I,fj^)  for  constant  k. 
Let  ( 1-v)  *1^  =  IjL 
From  equation  (11) 

(1+v)  *  I2  =  (l-Vo)l2 


Using  equation  (23)>  we  obtain 


(d)  = 


2(li)^+F[l-v,(l2)^] 


G(l-v,t^) 


(12)^-  (1-Vq)  G(l+v,t^; 


(26) 


(27) 


where  the  subscript  n  designates  the  values  at  ^  =  t^ . 
Hence,  l2(t)  can  be  evaluated  numerically.  Let 

Ax  =  —  ,  where  m  is  the  total  number  of  equal  space 

m  r- 

steps.  Define 

=  (l-l)Ax 

Pj  =  p(tj) 
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The  first  and  second  partial  derivatives  of  0  with 
respect  to  x  can  be  approximated  by  the  central  differ¬ 
ence  formulae: 


^ij  "  sfec  (  ®i-l,j  ®i+l,j) 


®ij  "  ^®i-l,j"  ^®ij^  ®i+l,J^ 


Equations  (12)-(15)  can  be  written  as  the  following 
matrix  difference  equations; 
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in 


t  ) 


(40) 


In  our  numerical  analysis,  one  extra  space  station 
at  X  =  A  +  Ax  is  added  for  the  convenience  of  repre¬ 
sentation.  Equations  (39)  and  (40)  contain  terms  nonlinear 
in  0^^  and  An  interative  method  is  used  to  solve 

this  non-linear  problem,  i.e.  we  assume  the  values  for 
0^^  and  and  evaluate  0^^  and  by  taking  s^^  and 

T,  as  known  functions.  Such  iterations  continue  until 
in 

a  certain  convergence  criterion  is  satisfied.  If  is 

known,  we  evaluate  y^^  by  the  numerical  technique  given 
in  [7] .  Let 

^in  "  ^in  ^i4-l,n  '^in  ^ 


where  is  a  two-by-two  matrix  and 

vector . 

From  equation  (28),  we  have 
^^in^i-l,n'^®in^  ^in  ^  *^in  ^i+l,n  “  ^in” 


is  a  column 
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Therefore, 
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“^in  "  ^in^^in"^in'^i-l,n^ 


where 


^in  “  ^^in^i-l,n  ®in^ 
Since  y^^  =  0 
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Pm  =  ^in  =  ° 


(^^5) 


Equations  (42)  and  (43)  are  the  recurrence  formulae  for 
Pfn  and  Therefore,  the  value  of  and  can 

be  determined  for  1  =  1,2, ...m+1.  Substituting  equation 
(4l)  Into  equation  (29)  we  can  show  that 

ym+2,n  ■  (Vm+l,n«n)'4H„-EAn-RnVl,n)  (“S) 


where 


R  =  E  P  +  P 
n  n  mn  n 


(47) 


After  yjjj^.2  n  determined  by  equations  (46),  y^^(  1=1,  2,  .  .  .m+1 ) 
and  p  can  be  evaluated  by  equations  (4l)  and  (l6)  respectively. 


In  the  case  where  p  Is  prescribed  as  a  constant.  It  is 
found  by  the  following  procedure  that  the  external  pressure 
relaxes.  Multiplying  both  sides  of  equation  (12)  by  dx 
and  integrating  from  x  =  0  to  x  =  X,  we  have 

A 

A^Ig  *  9'(A)  =  -  2  +  J  x$(9-x)dx  (48) 

o 


o 

Again,  if  we  multiply  both  sides  of  equation  (12)  by  x  dx 
and  Integrate,  we  obtain 


*  [A^9'(A)  +  8p] 
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-  ^  A^  +  -^  r  x^<I>(9-x)  dx 

^4 


(^9) 
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Eliminating  *  &'(^) 


from  equations  (48)  and  (49)  and 


solving  for 


P  = 


p,  we  have 
X 

[8pl2  +  y  x(l+  |)(l-  ^)  <I'(e-x)dx] 
o 


(50) 


hence,  the  relaxing  external  pressure  at  any  time  can  be 
determined . 

In  our  numerical  calculation,  we  computed  I^(t), 

l2(t),  ''(4)  by  equations  (22),  (27),  (21)  and 

(9)  respectively.  These  values  are  plotted  as  functions 
of  t  in  Figure  3.  The  functions  0  and  <i‘  were 
evaluated  by  iterative  procedure  from  equations  (39) >  (40) 
(4l),  (46),  etc.  If  p  was  prescribed,  the  calculation 
continued  until  the  absolute  value  of  the  difference  between 
the  value  of  p  and  the  average  value  of  p's  found  in 
the  five  previous  iterative  processes  was  less  than  0.001. 
If  p  was  prescribed,  the  convergence  criterion  used  here 
was  that  the  absolute  value  of  the  difference  between  the 
value  of  p  calculated  by  equation  (50)  and  the  average 
value  of  p's  found  in  the  five  previous  interative  pro¬ 
cesses  was  less  than  0.001.  In  our  computation,  we  used 
a  log  time  step,  therefore,  the  variations  of  Q  and 
at  very  early  times  could  be  taken  into  consideration.  The 
calculated  values  of  G  and  <I>  at  each  time  step  were 
stored  for  the  sake  of  the  history  dependence  of  the 
function  F  [equation  (24)].  For  the  purpose  of  comparison, 
numerical  calculations  were  carried  out  for  X  =  6  and 
p  =  0.8  (l.e.  80^  of  the  classical  buckling  pressure  of 
a  complete  sphere)  for  different  time  intervals  and  space 
intervals.  The  results  are  shown  in  Figure  4.  Although 
there  are  discrepancies  between  these  curves  for  the 
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different  space  intervals  used,  the  critical  times  obtained 
from  these  curves  show  no  noticeable  difference.  We,  there¬ 
fore,  used  Ax  =  0.125  and  A( log^Qt )  =0.5  in  the  subse¬ 
quent  computations. 

When  t  approaches  p  increases  very  rapidly. 

The  iterative  process  becomes  divergent  as  t  exceeds  a 
certain  value  This  divergence  of  our  iterative  pro¬ 

cedure  arises  from  the  fact  that  there  does  not  exist  a 
equilibrium  state  in  the  neighbourhood  of  the  equilibrium 
state  Just  before  buckling  defined  by  positive  Increaments 

of  time.  We  define  t  such  that  at  t  =  t  ,  the 

cr  cr 

logl0tcr+  0-005 

Iterative  process  converges  and  at  t  =  10  = 

1.012  tgjpj  the  iterative  process  diverges.  The  maximum 

relative  error  of  t  found  from  this  definition  is  less 

cr 

than  1.2^. 

Numerical  Results  and  Discussions 

A  program  for  numerical  computations  was  written  in  the 
BALGOL  machine  language  and  all  calculations  were  made  on  the 
IBM  7090  digital  computer  at  Stanford  University  Computation 
Center.  These  facilities  were  in  part  made  available  under 
the  National  Science  Foundation  Grant  NSF-GP  9^8.  The 
average  vertical  deflection  under  different  pressures  are 
plotted  against  time  for  A  =  6  and  A  =  7  in  Figures 
5  and  6.  The  critical  times  under  different  pressures  are 
shown  in  Table  1  and  Figure  7. 

From  Figure  'J ,  we  find  that  the  critical  time  decreases 
when  the  pressure  increases.  The  rate  of  decrease  becomes 
very  rapid  when  the  pressure  is  close  to  the  elastic  buckling 
pressure.  As  we  know,  the  critical  time  is  zero  when  the 
external  pressure  reaches  the  elastic  buckling  pressure  based 
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on  the  mod'.'li  for  instantaneous  loading.  Since  the  elastic 
buckling  pressure  for  A  =  7  is  higher  than  that  for 

A  =6^^^,  under  the  same  external  pressure,  the  shell  with 
A  =  7  possesses  longer  critical  time  than  the  shell  with 
A  =  6.  Figures  8,  9  and  10  show  the  distributions  of  B, 
and  at  different  times  for  A  =  6  and  p  =  O.85. 

The  Q  curves  in  Figure  8  are  the  slope  distributions  be¬ 
fore  buckling.  After  buckling,  the  values  of  9  are  pre¬ 
sumably  positive  in  the  whole  range  of  x.  Before  buckling, 
the  vertical  deflection  grows  with  time  but  the  shape  of 
the  deflection  curve  remains  approximately  the  same.  From 
Figures  9  and  10,  we  find  that  the  moments  near  the  shell 
edge  relax  and  that  the  moments  in  the  region  outside  the 
boundary  layer  grow  with  time  as  a  result  of  magnification 
of  the  slope  curve  (Figure  8) . 

The  case  in  which  p  =  I.O92  was  evaluated  and  the 
corresponding  relaxation  curve  of  external  pressure  p  is 
shown  in  Figure  11  for  A  =  6.  Figures  12,  13  and  l4 
represent  the  distribution  of  9,  m^^  and  m^  at  different 

time  stages.  When  time  Increases,  the  deflections  in  the 
central  portion  of  the  shell  decrease  but  the  deflections 
elsewhere  Increase  so  the  average  vertical  deflection 
remains  a  constant.  In  Figure  13  and  14,  moments  decrease 
with  increasing  time  because  of  relaxation  of  external 
pressure . 
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logio' 

[t  ) 

^  cr^ 

p 

A  =  6 

II 

0.80 

-0.22 

0. 14 

0.85 

-0.59 

-0.17 

0,90 

-1.05 

-0.48 

0.95 

-1.75 

-0.91 

Table  1 

Critical  Times  for  A  =  6  and  A  =  7 
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FIG.  1. 


Geometry  of  Clamped  Shallow  Spherical  Shells 


PIQ.  2.  Average  Vertical  Deflection  of  Clamped 
Shallow  Spherical  Shelia 
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FIG.  4.  Comparison  of  Average  Vertical  Deflections  for 
Different  Time  Intervals  and  Different  Space 
Intervals 
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FIG.  t.  Average  Vertical  Deflections  for  X  «  7 
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PIG.  11.  Relaxation  of  External  Pressure  When  p  Is  Specified. 


FIG.  13«  >"1.  Distributions  When  p  Is  Specified 


